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SECTION  1 


INTRODUCTION 


Bird  impact  is  a  critical  problem  in  the  latest  production  United 
States  Air  Force  aircraft,  the  F-I6  and  the  F/FB-111,  for  its  occurrence 
is  not  infrequent:  415  incidents  were  reported  in  the  period  1963-1972 
involving  the  loss  of  five  U.S.  planes.  In  the  F-16,  bird  strike 
protection  is  provided  by  allowing  the  canopy  to  deform  plastically  to 
absorb  the  impact  energy.  However,  there  is  a  possibility  that  the 
canopy  may  deform  sufficiently  to  result  in  contact  with  the  pilot's  helmet. 

It  was  the  purpose  of  this  study  to  develop  a  simulation  tool  which 
will  enable  the  designer  to  simulate  bird/canopy  impact  and  any  subse¬ 
quent  contact  of  the  canopy  with  the  pilot.  In  this  simulation  tool, 
the  pilot  is  represented  by  a  complete  head-spine  model  including  a 
helmet-head  interaction  model  and  a  skull-brain  model.  The  severity  of 
the  environment  and  the  degree  of  pilot  impairment  or  injury  can  be 
evaluated. 

For  the  purpose  of  simulating  the  pilot,  the  simplified  spine  model 
described  in  Belytschko  and  Privitzer  (1978)  AMRL-TR-78-7  was  used. 

In  order  to  permit  the  evaluation  of  bird  impact  in  terms  of  injury 
and  impairment  of  the  pilot,  a  skull-brain  model  based  on  the  maximum 
strain  criterion  (MSC)  of  McElhaney,  et  al  (1973)  was  incorporated  in 
the  model.  These  additional  developments,  the  helmet  model,  and  its 
interaction  with  the  head  are  described  in  Section  3. 

The  basic  simulation  tool  chosen  for  the  canopy  is  a  large 
deflection  nonlinear  finite  element  code  for  the  transient  analysis  of 
plates  described  by  Marchertas  and  Belytschko  (1974),  and  Belytschko 
and  Marchertas  (1974)  which  was  incorporated  in  the  general  purpose, 


program  Structural  Analysis  of  Man  (SAM) .  A  contact-impact  algorithm 
was  developed  for  modelling  bird-canopy  impact  and  the  impact  of  the 
canopy  with  the  pilot's  helmet.  The  resulting  simulation  tool  is 
described  in  Section  2. 

A  third  development  necessary  for  this  project  was  the  coding  of 
appropriate  material  models  for  the  canopy.  The  canopy  material  is 
nonlinear  but  it  exhibits  little  permanent  deformation.  Therefore,  a 
bilinear  elastic  law,  which  reproduces  the  initial  stress  strain  curve 
of  the  polyurethane  material,  but  springs  back  to  its  initial  shape  was 
formulated.  This  material  law  and  its  fit  to  experimental  data  are 
described  in  Section  4. 

Section  5  gives  the  results  of  our  simulations  of  bird  strike.  The 
initial  simulations  were  performed  without  a  pilot  in  order  to  study 
the  effectiveness  of  the  model  in  reproducing  experimental  results  for 
canopy  deformations.  It  was  found  that  only  a  coarse  mesh  could  be 
used  in  a  complete  simulation  if  core  storage  and  computer  cost  limitations 
were  to  be  met,  but  this  mesh  was  only  marginally  satisfactory  in 
reproducing  observed  behavior  of  the  canopy.  One  simulation  of  the 
complete  event  was  made  to  demonstrate  its  feasibility,  but  in  the 
remaining  simulations  which  involved  the  pilot,  the  bird  impact  was 
represented  simply  by  an  impulse  at  the  strike  point  with  a  coarse  mesh. 

A  computer  graphics  capability  was  developed  to  display  the  response  of 
the  canopy  and  pilot,  and  it  proved  quite  valuable  in  assessing  the  results. 

The  simulations  show  that  the  injury  potential  is  quite  sensitive 
to  the  initial  impulse,  its  locations  relative  to  the  pilot,  and  the 
pilot's  location  relative  to  the  canopy.  One  simulation  with  a  8  cm 
maximum  deflection  of  the  canopy  yielded  a  maximum  strain  criterion 


(MSC)  which  was  certainly  tolerable.  On  the  other  hand,  a  33%  increase 
in  the  impulse  yielded  MSC  values  associated  with  irreversible  injury. 
These  studies  were  only  preliminary,  in  that  certain  constants,  such  as 
the  helmet-head  spring  constant  values  were  not  available  while  most 
of  the  studies  were  conducted.  However,  the  results  indicate  that  this 
program,  after  validation,  should  provide  a  valuable  tool  for  the 
study  of  bird-canopy  impact.  Furthermore,  the  addition  of  plate  element 
to  the  SAM  model  should  be  useful  in  the  simulation  of  crash  environment 
where  they  can  be  used  to  represent  surfaces  which  impact  with  the  pilot 


SECTION  2 


COMPUTATIONAL  PROCEDURE 
2 . 1  Definitions  and  Coordinates 

The  mathematical  model  represents  the  structure  by  a  collection  of  rigid 
bodies  and  deformable  elements  in  a  three  dimensional  space.  The  format  of 
matrix  structural  analysis  and  finite  elements  is  used,  so  that  the  topology 
of  the  system  is  quite  arbitrary,  and  complex  structures  may  be  modeled. 

The  deformable  elements  are  interconnected  through  nodes,  as  shown 
in  Fig.  2.1. 

Two  types  of  nodes  are  used: 

1)  primary  nodes,  each  of  which  has  six  degrees  of  freedom  consisting 
of  three  translations  and  three  rotations;  located  at  the  centroids  of  the 
rigid  bodies; 

2)  secondary  nodes,  each  of  which  is  connected  through  a  rigid  body 
to  a  primary  node;  a  secondary  node  has  no  independent  degrees  of  freedom. 

An  arbitrary  number  of  secondary  nodes  may  be  associated  with  any  rigid 
body,  and  they  are  used  to  connect  deformable  elements  to  a  rigid  body  at 
points  other  than  the  centroid. 

In  order  to  describe  the  system,  we  will  define  three  types  of 
coordinate  systems: 

1)  a  fixed,  global  set  of  coordinates  (x,y,z),  or  x^; 

2)  body  coordinates  (x,  y,  z)^.;  a  triad  of  unit  vectors  b^,  $  and 

“►  _  _  _ 

b31  are  associated  with  each  node  I,  and  x.  y.  and  z  coincide  with  b^,  b9I 
and  ^3j>  respectively.  The  origin  of  the  x^  coordinate  system  must  be 
at  the  center  of  mass. 

3)  element  coordinates  (x,  y,  z);  a  set  of  element  coordinates  is 
associated  with  each  element,  and  the  element  coordinates  rotate  and  translate 
with  the  element  in  a  manner  to  be  specified  later.  The  x,  y,  and  z  axes  are 
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orientations  of  the  primary  nodes.  The  original  position  of  node  I  is  denoted 


by  (i=l  to  3  representing  the  x,  y,  and  z  components);  the  new  positions 


x^  are  obtained  by  adding  the  displacements  to  the  original  position,  so 


Xil  =  Xil  +  Uil 


(2.6) 


The  velocities  and  accelerations  of  the  primary  nodes  are  denoted  by  and 
respectively. 

The  orientation  of  each  primary  node  is  described  by  a  triad  of  orthogonal 

"►  -f 

unit  vectors  b^,  °2j_’  b3i*  whicb  rotates  with  the  node.  The  angular  velocities 
and  angular  accelerations  are  specified  in  body  coordinate  components, 
and  respectively. 

The  inertial  properties  of  the  primary  nodes  are  described  by  the  masses  of 
primary  nodes,  p^.,  and  the  principal  moments  of  inertia  of  the  primary  nodes, 

III ,  *31'  *be  n0<*a*  vectors  b^  must  coincide  with  the  principal  coordinates 

of  the  moment  of  inertia  tensor. 

The  forces  and  moments  at  the  nodes  are  denoted  by  and  M^,  respectively, 

£XC 

and  may  be  subdivided  into  externally  applied  forces  and  moments,  and 

—  6Xt 

M  ,  and  the  forces  and  moments  due  to  the  resistance  of  the  deformable 
elements,  F^*1  and  Moments  are  always  treated  in  the  nodal  coordinate 

systems,  x^. 

The  element  quantities  are  extracted  from  the  global  quantities  in  the 
usual  manner  of  maxtrix  structural  analysis  by  a  Boolean  connectivity  matrix 

(p) 

,  so  that 


j(e) 

iA 


u 

AI  il 


(2.7) 


where 


(e) 

A;  =  1  if  the  Ath  generic  code  of  element  e  corresponds  to  the 

Al 


*(e) 

lAI 


Ith  primary  node  of  the  system. 
0  otherwise 


The  contributions  of  the  element  nodal  forces,  and  {n/6^} 

nodal  forces  are  obtained  by  the  equation  conjugate  to  Eq.  (27) 


to  the  total 


int 

=  y 

»(e) ,(e) 

il 

L 

e 

AI 

iA 

int 

=  Y 

il 

L 

A I 

iA 

(2.8a) 

(2.8b) 


e 

Since  the  connectivity  natrix  is  Boolean,  the  operations  of  Eqs.  (2.7)  and 
(2.8)  are  never  carried  out  as  matrix  multiplications.  Equation  (2.7)  simply 
represents  the  extraction  of  a  smaller  array  from  larger  array,  while  Eq.  (2.8) 

represents  the  addition  of  a  smaller  array  into  the  appropriate  slots  of  the 
larger  array. 


2 . 2  Equations  of  Motion 

The  equations  of  motion  for  the  system  are  written  separately  for  the 
translational  and  rotational  degrees  of  freedom.  The  three  translational 
equations  are 

PIU1I  =  F±I  -  Fix  (no  sum  on  I)  (2.9) 

Since  the  x^  coordinates  are  principal  coordinates,  the  rotational  equations  are 


t  I  /  T  ,  .  r^Gxt  ,,int 

11T  a  +  (I„T  -  1„T)  co  co  =  M  _  -  M  _ 

II  xl  v  31  21  yl  zl  xl  xl 

T  _L/T  T  \  rGXt  ^int 

I0T  a  +  (I  -  I  )  co  co  =  M  -  M  _ 
21  yl  II  31  xl  zl  yl  yl 


t  —  .  t  n  —  —  wext  rant 

I,T  a  +  (I.  -  I  )  co  co  =  M  -  M  _ 

31  zl  21  II  xl  yl  zl  zl 


(2.10) 


(no  sum  on  I) 


On  an  element  basis,  the  nodal  forces  and  moments  are  first  computed 
in  terms  of  the  corotational  element  coordinate  system  x^.  If  an  element  is 
connected  to  a  primary  node,  the  nodal  forces  are  then  transformed  into  the 
global  system  by 

(f(e))  =  [e]T  {f(e)}  (2.11) 

while  the  moments  are  transformed  into  the  nodal  coordinate  system  by 
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{m(e)}  =  [b]  [e]T{m(e) }  (2.12) 

The  resulting  nodal  forces  are  then  summed  to  give  the  total  Internal  forces 


and  moments. 

If  an  element  is  connected  to  a  secondary  node  A  which  is  associated  with 
a  primary  node  I,  then  it  is  necessary  to  account  for  the  rigid  linkage  in 
transforming  the  element  nodal  forces.  We  note  that  the  nodal  velocities  of 


the  secondary  node  are  given  by 


where 


where  double  subscripts  designate  a  difference  in  node  coordinates, 


(2.13) 


(2.14) 


e.  g. 


(2.15) 


It  can  then  be  shown  (Belytschko,  et  al.,  (1976))  that  the  nodal  forces  at 
the  primary  node  are  given  by  Eq.  (2.11),  while  the  nodal  moments  are  given  by 

{m(e)]  =  [b]  [e]T(i(e)}t+  [fi]T[b]  [e]T  {f(e)}  (2.16) 

Thus  while  for  primary  nodes  the  nodal  moments  are  simply  related  by  the 
coordinate  transformations,  for  secondary  nodes,  an  additional  moment  is 


introduced  in  the  transformation  because  of  the  moment  arm  between  the  nodes. 


2 . 3  Beam  Element 

The  beam  element  is  developed  in  terms  of  a  corotational  formulation  as 
described  by  Belytschko,  et  al  (1978).  In  a  corotational  formulation, 
displacements  and  rigid  body  rotations  of  each  element  may  be  arbitrarily 
large,  but  the  deformation  of  each  element  relative  to  the  corotational 
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coordinates  must  be  small.  This  does  not  limit  the  formulation  to  small 
deformations  but  imposes  the  requirement  that  the  number  of  elements  used  in 
the  model  be  sufficiently  large  so  that  the  deformation  of  each  element 
relative  to  its  rigid,  Cartesian  corotational  system  be  small. 

Standard  engineering  beam  theory,  often  known  as  Euler-Bernoulli  beam 
theory,  will  be  used.  This  theory  assumes  that  plane  sections  normal  to 
the  midline  remain  plane  and  normal.  As  a  consequence,  the  midline  of 
the  beam  must  rotate  with  the  nodes;  if  two  elements  are  connected  in  a 
straight  line,  the  midline  must  remain  continuous  and  continuously  differ¬ 
entiable.  Warping  of  the  beam  cross-section  will  be  neglected. 

Consider  a  generic  beam  element  with  nodes  1  and  2  as  shown  in  Figure  2.2. 

AAA 

The  corotational  coordinate  system  (x,  y,  z)  for  the  element  is  defined 

A 

by  requiring  the  x  axis  to  always  pass  through  the  two  nodes  with  its 

A 

origin  at  node  1.  The  y  axis  is  defined  by  the  average  rigid  body  rotations 

A  A  A  A 

about  the  x  axis,  and  the  z  axes  is  defined  by  its  normality  to  the  x-y 


plane. 

^dc f 

The  deformation  displacements  are  then  u^  ’  Uy  ’  Uz  ’  anc^  t*ie 

Ade  f 

rotation  0^  .  From  the  definition  of  the  corotational  coordinate  system, 

a  dcf  Ade f  Adef 

it  follows  that  u  vanishes  at  node  1,  while  u  and  u  vanish  at 
x  y  z 

both  nodes. 

The  deformation  displacement  fields  are  chosen  to  be  polynomials  and 


the  lowest  order  polynomial  consistent  with  Euler-Bernoulli  beam  theory 

A(jg  f  *ds  f 

is  used.  Since  the  transverse  displacements,  u  and  u  ,  must  be 

y  z 

continuously  differentiable,  cubic  polynomials  are  used  for  these  defor- 

AJgf  Mgf 

mation  displacements.  The  axial  displacement,  u  and  the  rotation  0 


IM> 


UNDEFORMED 


Figure  2.2  Beam  Element 
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’7’ 


X  > 


need  only  be  continuous,  so  a  linear  field  is  used.  Thus 


Adef  Adef  _ 
Ux  =  Ux2  5 


“  6J 

.2  .  _ 3 


(s  -  2?  +  c3)  Szl  + 1  <-r  +  ?3> 


3  A 

'  e 


z2 


Adef 
u  = 

y 

Adef  2  3  A  2  3  A 

u“et  -£(-?  +  2?z  -  eyl  +  «  (5  -  ?J)  ey2 


(2.17) 


9def  =  e  ?1  ? 

x  x21 


5  =  x/£ 

where  l  is  the  length  of  the  element  and  0^^  the  relative  rotation 

between  nodes  2  and  1.  The  rotations  9  ,  and  9  are  defined  by 

yl  zl  J 


zl 


yi 


Adef  .a 
du  (xT) 
y  I 

A 

dx 


_  a  def 
du 

_  z _ 

A 

dx 


I  =  1,2 


The  deformation  nodal  displacements  are  then 
T 

1  A  A  A  A  A 

Adef  \  =  6,0,.  9  .  ,  0  .  ,  0  ,  9  , 

-d  J  L  x21,  yl  zl’  yz  z2. 


(2.18) 


(2.19a) 


6  =  l  ^Ux21  X21  +  Uy21  y21  +  Uz21  Z21  +  2(ux„.  +  u2  +  u2  3  ^  (2>19M 

zi  yzi  zzi 


X21=  x2  *  X1 


ux21  =  ux2  ‘  Uxl 


(2.19c) 

(2 . 19d) 


Ade  f 

where  6»  u^  is  the  elongation  of  the  beam. 


A  A  A 

For  the  purpose  of  computing  0^^,  eyI>  and  ezI  in  the  deformed 


->o 


configuration,  nodal  components  of  e.^  =  and  e =  s^  in  the  undeformed 
configuration  are  stored  for  both  nodes  of  the  element.  The  nodal  (body) 


components  of  the  vectors  r^  and  s^  are  considered  invariant,  so  these 
vectors  rotate  with  the  nodes.  Consequently,  the  angle  between  r  and  e^  in 
the  deformed  configuration  gives  the  magnitude  of  rotation  of  the  midline 
relative  to  the  corotational  x-axis  of  the  element;  if  r  and  e^  are 

A 

coincident  in  the  deformed  configuration,  the  deformation  rotations,  d  and 

A  A 

9  vanish  for  the  element.  For  small  deformations,  the  magnitudes  of  0  and 
z  y 

A 

may  be  found  by 


/S  A 

9  e„  +  0 


■  -  e  =  e.  x  r 
y  2  z  3  1 


(2.20) 


The  element  components  of  r  may  be  found  by  Eqs.  (2.3)  and  (2.4)  and 
expressing  the  right  hand  side  in  element  coordinates  we  find 


$  e  +  ^  e  =  det 
y  2  z  3 


2  3 

0  0 

r  -r 


-r  +  r  e„ 
z  2  y  3 


(2.21) 


Thus,  by  equating  components,  it  follows 


A  A 

9  **  -r 

y  2 


A  A 

9  =*  r 

z  V 


(2.22) 


The  torsional  deformation  rotation  is  found  by  taking  the  cross-product  of 

— >  — >  A 

s^  with  Sj,  and  projecting  the  result  on  the  unit  vector  along  the  x-axis. 


which  gives 


x21 


“■>  , 

ej  •  x  s2)  -  det 


s,  s„  -  s„  s, 
ly  2z  2y  lz 


1 

CM 
f  Q> 

-> 

e3 

lx 

siy 

*lz 

2x 

S2y 

S2z 

(2.23) 


17 


The  strains  in  a  corotational  formulation  are  given  by  (Belytschko 


and  Hsieh  (1973)) 


Adef 

hu 

X 

A 

hx 


(2.24) 


For  an  Euler-Bemoulli  beam,  the  assumption  that  normals  to  the  midline 
remain  straight  and  normal  gives 


Ade  f 
u 

x 


AAA 

(x,  y,  z) 


Adef  *  „  A  &u2 

ux  (x,  0,  0)  -  y  — 

hx 


def 


A  def 
Ou 

_ 1 _ 

A 

&X 


Substituting  the  above  into  Eq .  (2.24)  gives 

2  Adef  . A 


Ade  f  a 

5u“er  (x,0,0)  A 

~  'a - -  y 

hx 


b  U  (x,0,0) 
z 


bx 


pAdef  A 

a  6  u  (x,  0,0) 

.  Z_^L_ - 

6x 


(2.25) 


For  the  purpose  of  determining  the  shear  strains  due  to  torsion,  we  consider 

A  A 

polar  coordinates  r,  0,  in  the  y,  z  plane.  Then 


2e 


r9 


*def 

hu 

r 


.  Adef 

_^e _ 

hr 


(2.26) 


Ade  f  Ade  f 

We  assume  u^_  =0  and  that  u^  varies  linearly  with  r;  this  assumption  is 

consistent  only  with  cylindrical  cross  sections.  Equations  (2.17)  and  (2.26) 
and  the  above  assumption,  along  with  our  neglect  of  warping,  then  gives 


re 


x21 


(2.27) 


r0  l 

Substituting  Eq.  (2.17)  into  Eq.  (2.25),  we  then  obtain  the  strain  displace¬ 
ment  equation  for  the  element  to  be 


‘re 


7  f1  0 


z(6?  -4)  y(4  -  6?)  z(6?  -  21  y(2 


L°  r 


;  i  H 


(2.28) 


18 


The  above  relation  defines  the  strain-displacement  matrix  [b],  Using  the 
theory  of  Belytschko  and  Hsieh  (1973),  it  follows  then  that  the  deformation 
nodal  forces  are  given  by 


l  M 


x  i  dv 


(2.29) 


when  V  is  the  volume  of  the  element  and  where  the  deformation  nodal  forces 
are  conjugate  to  the  deformation  nodal  displacements  in  the  sense  that  their 
scalar  product  gives  work,  so  that 


,  -def  *  a  a  a  a 

x2,  x2,  yl,  zl,  y2,  z2J 


(2.30) 


In  order  to  minimize  computational  expense,  it  is  often  desirable  to 
circumvent  the  volume  integration  of  Eq.  (2.29)  by  defining  generalized 
stresses 


J.  r  5 


a  „  dA 


(2.31a) 


i  z  a  dA 
JA  x 

r  ,  „ 

-  i  y  a  dA 

"A  x 


(2.31b) 


(2.31c) 


where  A  is  the  cross-sectional  area  of  the  beam.  The  conjugate  generalized 


strains  are 


^2Adef 

9  u 


(2.32a) 


^2  Adef 
a  u 

z 


(2.32b) 
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For  a  linear  material,  these  definitions  of  generalized  stresses  and  strains 
may  be  combined  with  Eqs.  (2,17),  (2.29)  and  (2.31)  to  yield 


' 

m 

yi 

i--CKb]  ■ 

r 

A  1 

V 

A 

\  .  _  r  b-,  \ 

1  K  1 

r  •  ~^j 

A  i 

Vi 

m 

L  y2. 

1  1 

y2  J 

y  I 

1 

r 

CD> 

,/U 

ir"z2! 


where 


[Kbl  »  kb 

l  l 


El. 


f- 

9Z1 

_ez2. 

[Kb]  < 

2  ! 

1 

i 

2  -a>. 

i 

). 

l 

4+0^ 

zl  V 


z2 


where  E  is  Young's  modulus,  I  are  the  section  moduli, 


(2.33a) 


(2.33b) 


(2.34a) 


(2.34b) 


I 

yy 


r  a2 

z  dA 
‘  A 


I 

zz 


f  2 

=  y  dA 
'A 


(2.35) 


and  c^  is  a  stiffness  proportional  damping  factor.  The  shear  factor  is 
given  by 


24  (H-u)  I 

A  l2 

s 


(2.36) 


where  v  is  Poisson's  ratio  and  A  is  the  effective  area  in  shear. 


s 

While  the  bending  stiffness  and  torsional  stiffnesses  have  here  been 
expressed  in  terms  of  elastic  material  constants  for  homogeneous  isotropic, 
elastic  beam,  in  most  elements  of  the  biodynamic  model,  the  cross  section  of  the 
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elements  is  quite  inhomogeneous  and  anisotropic,  and  the  overall  stiff- 
3  b  t 

nesses  k  ,  k  and  k  ;  instead  the  material  moduli  have  been  determined 
experimentally.  Therefore,  these  overall  stiffnesses  are  used  directlv  in 
the  model. 

Once  the  nodal  deformation  forces  are  known,  the  nodal  forces  required 
for  the  equations  of  motion  can  be  obtained  in  two  steps.  First,  the 
remaining  six  nodal  forces  (the  beam  element  has  two  nodes,  each  with  6 
degrees  of  freedom)  are  found  by  equilibrium  requirements 


A  A 

f  .  =  -  f  0 
xl  x2 


A 

.  +  m  - 
zl  z2 


yi 


zl 


A  A 

m  ,  +  m  „ 

.  ..yi,  y.2 


A  A 


m  1  = 
xl 

-mx2 

A 

A 

f  = 

-  f 

y2 

yi 

A 

A 

f  = 

-  f 

z2 

Zl 

(2.36) 


The  global  components  are  then  found  by  Eq.  (2.3)  and  (2.4).  If  the  node 
is  a  slave  node,  Eqs.  (2.11-13)  are  used. 


2.4  Plate  Element 

Initially,  we  considered  solving  the  canopy-pilot  interaction  problem 
by  using  the  programs  for  the  pilot,  SAM,  and  the  plate  finite  element 
program  developed  at  Argonne  by  Marchertas  and  Belytschko  (1974),  SADCAT, 
interactively.  However,  we  found  that  considerable  core  storage  would  be 
required  to  keep  both  programs  in  core,  while  an  overlay  would  involve  too 
much  computer  time,  since  the  complete  simulation  involves  several  hundred  time 
steps,  and  the  two  programs  would  have  to  be  interchanged  in  fast  core  for  each 
time  step.  Both  SADCAT  and  SAM  were  written  in  terms  of  the  same  nodal  updating 
format  and  logical  structure,  so  it  was  decided  that  the  best  course  would 


be  to  simply  incorporate  the  plate  element  as  another  element  in  the  SAM 
library. 

The  plate  element  is  a  flat,  triangular  element  as  shown  in  Fig.  2.3. 
Although  a  curved  shell  such  as  a  canopy  must  be  approximated  as  piecewise 
flat  by  these  elements,  past  experience  (Marchertas  and  Belytschko  (1974)) 
indicates  that  it  effectively  solves  curved  shell  problems.  The  plate 
element  also  uses  the  corotational  formulation.  The  (x,  y)  plane  is  defined 
to  always  coincide  with  the  plane  defined  by  the  3  nodes  of  the  element;  the 
origin  of  the  (x,  y)  coordinate  system  coincides  with  node  1,  and  the 
x-axis  is  aligned  so  that  it  bisects  the  angle  between  sides  12  and  13;  see 
Fig.  2.2. 

The  plate  element  also  uses  the  corotational  formulation.  The  (x,  y) 
plane  is  defined  to  always  coincide  with  the  plane  defined  by  the  3  nodes  of 
the  element;  the  origin  of  the  (x,  y)  coordinate  system  coincides  with 

A 

node  1,  and  the  x-axis  is  aligned  so  that  it  bisects  the  angle  between 
sides  12  and  13;  see  Fig.  2.3. 

The  deformation  displacement  field  is  assumed  to  be  linear  in  the  plane 
of  the  element  and  cubic  normal  to  the  plane,  and  is  expressed  in  terms  of 
the  nodal  deformation  displacements  by 

''l  r 

!  do  Co] 

>  ■  ;  do  [o] 

;  [o]  r»fi 

J 

where  [bm]  are  the  linear  triangular  shape  functions  for  the  membrane  (in¬ 
plane)  displacements  (see  for  example  Zienkiewicz  (1978))  and  [Nf]  are  the 


•  ddef 
_ 


22 


cubic  flexural  shape  functions  of  Bazeley,  et  al  (1965).  The  nodal 
deformation  displacements  are 


;def 


[«]  [o] 


<510  601  9,  e.  e„  e„  e„  e. 

12,  23,  31,  lx,  ly,  2x,  2y,  3x,  3y 


(2.38) 


where  6^  is  the  elongation  of  the  side  defined  by  nodes  I  and  J,  which 
is  given  by  Eq.  (2.19b). 

In  order  to  compute  the  nodal  rotations,  it  is  necessary  to  store  the 
body  (nodal)  components  of  the  vector  normal  to  the  plate  n  =  (nx,  n^,  n^) 
for  each  node  of  the  element.  Invoking  an  argument  identical  to  that  used 
to  obtain  Eqs.  (2.20-2.22),  we  can  write  that  for  small  9 


0x  el  +  ey  e2  =  e3  x  n 


(2.39) 


In  a  deformed  configuration,  n  is  expressed  in  element  coordinate-components, 
and  Eq.  (2.39)  can  then  be  shown  to  yield 


0  =  -n 

x  y 


0  =  n 

y  X 


(2.40) 


The  strains  in  a  corotational  plate  formulation  are  given  by 


9u 


def 


e  = 
x 


9x 


-  z 


. 2"def 
9  u 

z 

-2 

3x 


(2.41a) 


£  = 
y 


„  "def 
9u 

_ 

3y 


-  z 


2"def 
9  u 

z 

~  "2 

9y 


(2.41b) 


9u 


def 


2e 


xy 


3x 


9udef  32Gdef 

+  — - - 2z  - - - 


9y 


3x  9y 


(2.41c) 


Inserting  Eqs.  (2.37)  into  Eq.  (2.41),  we  obtain 


{£}  =  Jb®J  {6}  -  z  [Bf]  {§} 


(2.42) 


where 


(e}T  =  [e  ,  e  ,  2e  1 
*-  x’  y’  xy-1 


(2.43) 


on  - 


3N  ®/3i 
xl 

3N  ^/ay 

3N  ^/3y  +  9N  “j/Bi 


(2.44) 


f  32N  f/3x2 

[Bf] 

32N1f/3y2 

2  f  a  A 

3  N  /3x  3y 


(2.45) 


The  conjugate  nodal  deformation  forces  are  separated  similarly  into  membrane 
forces  and  moments.  The  membrane  forces  are  given  by 


[b®]t  {5}  dV 
*vT 


(2.46) 


where  (a)  is  the  corotational  stress  maxtrix 


{a}T  =  [a  ,  a  ,  a  1 
*-  x  y  xyJ 


(2.47) 


which  is  obtained  from  the  strains  by  the  stress-strain  law.  The  nodal 
deformation  moments  are  given  by 


{m}  *  "f  z  [Bf]T  {a} 

*'  XT 


(2.48) 
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where 


{m}  =  Kx’  mly*  m2x’  m2y’  m3x’  m3rl 

The  volume  integrals  of  Eqs.  (2.46-48)  are  evaluated  numerically  using  3 
points  in  the  (x,  y)  plane,  and  at  each  of  these  points  in  the  plane  a 
variable  number  of  points  in  the  z  direction  (through  the  depth  of  the 
plate  element)  specified  by  the  user  are  used  for  the  integration. 

The  transverse  nodal  forces  f  are  then  obtained  by  equilibrium 
requirements.  Moment  equilibrium  about  the  x  and  y  axes  gives 


-  — 

' 

L 

2z 

1 

*3  *3 

“lx  +  “2x  +  “3x 

'  =  2A 

f, 

3z 

-*2  -y2 

m.  +  +  m0 

L  ly  2y  3y  J 

and  equilibrium  in  the  z-direction  then  gives 


(2.49) 


(2.50) 


(2.51) 


where  A  is  the  area  of  the  element.  The  nodal  forces  are  then  transformed 
into  the  global  coordinates  and  utilized  in  the  equations  of  motion  as 
internal  nodal  forces. 


2 . 5  Time  Integration  and  Stability 

In  this  investigation,  time  integration  was  performed  by  the  Newmark 
8-method  with  8=0,  which  is  an  explicit  time  integration  method  with 
stability  and  convergence  properties  that  are  almost  identical  to  the 
central  difference  method.  The  Newmark  6-method  is  sometimes  more  con¬ 
venient  than  the  central  difference  method  because  all  nodal  variables 
and  their  time  derivatives,  i.e.  velocities  and  accelerations,  are  defined 
at  each  time  step,  whereas  the  central  difference  method  defines  velocities 
only  at  times  equidistant  between  time  steps,  often  called  half-time  steps. 

In  the  Newmark  8-formulas,  the  velocities  and  displacements  are 
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updated  in  each  time  step  by 


’  j+1  *  j  ,  1  ..  ...j  ,  ..j+1. 

Jil  =  uil  +  2  At  (uil+  uil  ) 


u^1  -  uj  +  At  u J  +  £  ut2  'u  .1 
ll  ll  ll  2  ll 


(2.52) 


where  At  is  the  time  step  and  superscript  dots  denote  the  time  step  number. 

For  the  rotational  degrees  of  freedom,  Eq.  (2.52)  cannot  be  used 
directly  because  the  orientations  are  described  by  the  nodal  (body) 
vectors  b^,  and  the  rates  of  these  vectors  are  not  equivalent  to  the 
angular  accelerations  and  velocities. 

The  counterpart  of  Eq.  (2.52)  for  the  unit  vector  b  is 


3 


*J+1 


-  j  *1 

biJ  +At“^ 


+  \  At2 


d2 


dt 


(2.53) 


d^J+1 

dt 


Si. 

dt 


♦J 


dt2 


-ir 


(2.54) 


From  vector  analysis,  it  follows  that 


db. 


dt 


U)  X  bJ 


(2.55) 


dS 


U)  x  (u)  x  bt)  +  (a  X  b^) 


(2.56) 


To  obtain  the  updated  components  of  the  unit  vectors,  we  temporarily  fix 

the  x^  coordinates  and  consider  a  particular  unit  vector;  we  then  dot 

Eq.  (2.53)  with  the  unit  vector  corresponding  to  the  component  of  b^  to 

be  updated,  and  substitute  Eqs.  (2.55)  and  (2.56)  into  the  resulting 

—  — > 

expressions.  For  example,  the  updated  x-component  of  b^is  found  by 
letting  i  *  3  in  Eq.  (2.53)  and  taking  the  scalar  product  with  b^,  which 
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yields  after  some  simplification  to 


b^+1  =  At  tu  ^  +  t  At2  (03  UJ  ^  +  a  j) 
3x  y  2  x  z  y 


Similarly 


b^+1  =  -At  03  j  A  t2  (03  ^  03  ^  -cr  J) 

3y  x  2  y  z  x 

b^+1  -  At  03  j  +  ^  &t2  (03  ^  03  ^  +  or  j) 
ly  z  2  x  y  z 


(2.57) 

(2.58a) 

(2.58b) 


Normality  and  orthogonality  of  the  unit  vectors  are  used  to  find  the  remain¬ 
ing  components.  It  is  assumed  that  the  rotations  during  a  time  step  are 
small  so  that  second  order  terms  can  be  neglected.  From  normality,  it 
follows  that 

2 

-  1 i, 

3y 


•  b  -  (Ht'Y  -  (^I1)2] 1/2 


(2.59) 


and  if  we  assume  that  b|^  ~  1,  orthogonality  yields 


p.+i  -  -(bJ+i  +  b \+i  bj+i)  bJ+i 

lz  \  3x  ly  3y  /  3z 
— i’t’l 

Finally,  the  component  b^x  can  be  found  by  the  normality  condition 


(2.60) 


* 1  ■  b  -  K'?  -  &:1)2 


1/2 


(2.61) 


The  vector  b£  is  then  found  by  the  cross-product. 

The  components  of  the  b^  and  b^  vectors  given  above  are  in  terms  of 
the  x^-coordinates  (nodal  coordinates)  at  time  step  j.  These  components 
are  then  transformed  to  global  components  by  the  transformation  (2.3). 

In  explicit  time  integration,  unless  the  time  step  is  sufficiently 


small,  numerical  instabilities  will  render  the  results  worthless.  Although 
stability  analyses  cannot  be  performed  for  the  fully  nonlinear  case. 


experience  indicates  that  a  linearized  analysis  gives  useful  estimates  of 
the  stable  time  step.  The  destabilizing  effects  of  nonlinearities  can 
usually  be  overcome  by  reducing  the  time  step  that  is  used  by  107.  to  20% 
below  the  linear  stability  limit.  Although  this  does  not  always  insure  a 
stable  computation,  stability  can  be  checked  after  the  computation  by 
examining  the  energy  balance. 

A  linearized  analysis  of  the  Newmark  B-method  indicates  that 

9  r  (  9'  1/2  -i 

itsu-  \  [l  +U-  )  -  pj  (2.62) 

max  L" 

where  u>  is  the  highest  frequency  of  the  system  and  p.  is  the  fraction 
max 

of  critical  damping  in  the  highest  frequency.  The  frequency  <o  corresponds 
to  the  maximum  eigenvalue  of  the  eigenvalue  problem 

[K]  jx}  -  ua2  [m  {x}  (2.63) 

where  [K]  is  a'  linear  global  stiffness  matrix  and  [Ml  is  the  mass  matrix. 

This  stability  limit  is  obtained  by  performing  a  Fourier  analysis  of  the 
uncoupled  equation 

•  2 

x  +  2  p.  u)  x  +  u>  x  =•  0  (2.64) 

No  provisions  are  made  in  the  program  for  directly  computing  the  maximum 
eigenvalue  w,  for  it  would  require  considerable  programming  and  computer 
time,  and  it  would  only  be  of  marginal  usefulness  in  nonlinear  problems 
since  [k]  and  [M]  can  change  dramatically  with  time.  Instead,  we  have 
developed  formulas  for  estimating  u>.  These  are  based  on  Rayleigh's  theorem, 
which  allows  us  to  bound  the  maximum  frequency  of  the  system  by  the 
maximum  frequency  among  all  of  the  elements  of  the  system  (see  Hughes, 


et  al  (1978)),  i.e 


uj  S  sup 
max 


uu 


ele 

max 


(2.65) 


0l0 

where  uj  ^  is  the  maximum  of  the  highest  frequency  of  all  elements  in 
the  system. 

The  maximum  frequency  of  any  of  the  elements  can  easily  be  obtained 
exactly  or  approximated  sufficiently  accurately  for  our  purposes.  First 

Q 

we  will  consider  a  spring  element  with  a  spring  constant  k  .  The  element 
stiffness  [k]  and  mass  [m]  matrices  are  given  by 


< - 1 

!  k 


— 

— 

1 

-1 

r  1 

m. 

0 

-1 

1 

,  mj- 

u 

1 

0 

m2 

(2.66) 


where  m^  and  m2  are  the  lumped  masses  at  the  nodes.  The  maximum  frequency 
is  given  by 


ou  -  -  (2.67) 

max  ml  +  m2 

Whenever  the  system  consists  of  a  system  of  springs  connecting  rigid 
bodies,  all  springs  in  parallel  between  any  two  rigid  bodies  should  be 
summed  to  obtain  the  total  spring  constant  between  the  two  rigid  bodies. 

For  the  purpose  of  using  Eq.  (2.67),  half  of  the  lumped  masses  at  each 
of  the  nodes  should  be  used.  While  in  principle  the  entire  lumped  mass 
could  be  used,  it  should  be  noted  that  the  same  mass  cannot  be  ascribed 
to  two  elements  in  estimating  the  maximum  frequency.  Therefore,  if  all  of 
the  mass  is  used  in  estimating  the  frequency  for  anyone  element,  the 
adjacent  elements  will  have  no  mass,  which  will  increase  the  frequency 
associated  with  that  element. 


When  the  lumped  masses  arise  from  the  density  of  the  element  itself, 


ml  “  ra2  =  2  P  A  ^ 


(2.681 


where  o  is  the  density.  If  the  stiffness  is  associated  with  an  elastic 
material  with  Young's  modulus  E,  then 


k  =*  AE/ 1 


Substituting  these  values  into  Eq.  (2.67),  we  find 


(2.69) 


2  4E 

U)  *  — =- 

max  H2 


(2.70) 


where  c,  given  by 


2  E 


(2.71) 


is  the  wavespeed  in  the  material.  Substituting  this  result  with  Eq.  (2.62), 
we  obtain 


tttS^  (  Vp,2  +  1  -y) 


(2.73) 


so  that  when  p  *  0,  the  stability  limit  on  the  time  step  corresponds  to 
the  time  for  the  elastic  wave  to  traverse  the  length  of  the  element. 

When  an  element  with  flexural  stiffness  connects  two  nodes,  its  effects 
on  stability  must  be  considered  separately.  For  this  purpose,  it  is 
sufficient  to  consider  bending  in  one  plane.  The  stiffness  matrix  (if 
the  shear  effect  is  zero)  is  then  (see  Przeminiecki  (1968)) 
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12  61 


[k]  =  ^ 

r 


-12  61 


- 61  U 


12  -61 


(2.74) 


symmetric 


and  we  will  consider  a  mass  maxtrix 


210  0  0  0 


210  0 


symmetric 


(2.75) 


The  maximum  frequency  is  given  as  a  function  of  a  by  the  roots  of  the 
equation 


o'  (p,  -  0. 115^)  +  o(0.229  -  8u)  +  12  =  0 

a  2 
pA  l  uj 

M.  _ 

420EI 

For  o  =  17.5,  which  is  the  lumping  we  use  whenever  the  mass  arises  from 
the  density  distribution  of  the  element  itself,  the  solution  is 


(2.76) 


192  c  l 
.„4 


(2.77) 


For  a  typical  model,  it  is  possible  that  either  the  frequency  associated 

with  axial  response,  Eq.  (2.67),  or  flexural  frequency,  Eq.  (2.77),  is 

the  largest.  In  choosing  the  time  step,  it  is  the  largest  frequency  which 

must  be  used.  In  comparing  Eq.  (2,71)  and  (2.77),  it  can  be  seen  that  the 

axial  frequencies  increase  inversely  with  l,  whereas  the  flexural  frequen- 

2 

cies  increase  inversely  as  i  .  Hence,  for  sufficiently  refined  models. 


the  flexural  frequency  will  govern  the  time  step. 
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The  frequency  associated  with  a  torsional  mode  is  given  by 


4kL 


max  i  +  i 

xl  x2 


(2.78) 


This  frequency  should  also  be  checked,  but  it  is  generally  smaller  than 
the  axial  frequencies. 

For  the  plate  element,  membrane  frequencies  can  be  estimated  by 
Eq.  (2.71)  with  the  smallest  distance  in  the  element  from  a  node  to  the 
opposite  side.  Equation  (2.77)  is  used  in  the  form 


2  2 

2  16  c  h 


(2.79) 


where  h  is  the  thickness  of  the  plate;  SL  is  used  in  the  same  sense  as  for 
the  membrane  frequencies. 

The  relation  between  stiffness  proportional  damping  c,  in  Eq.  (2.33) 


and  p  is  given  by 
c. 


y  = 


U) 

max 


(2.80) 


If  artificial  viscosity  is  used,  the  artificial  viscosity  constant 
corresponds  with  p.  The  formulas  for  stability  calculations  are  summarized 
in  Table  2.1. 

To  insure  post  facto  that  the  computations  are  stable,  energy 
balance  checks  were  made  as  follows:  The  kinetic  energy  at  each  time  step 
j  is  defined  by 

.  2 


Tj  =  2  pi  (u{i>  +  2  Iil  (wil* 


(2.81) 


where  I  is  summed  over  the  primary  nodes  in  the  system.  The  external  work 
is  defined  by  a  trapezoidal  integration  in  time,  so 
_,ext,o  _  . 


W 


(2.82) 


ext,j  ext,j-l  1^  .  j  _  (pext»J  + 

W  w  +2  tun  uiX  )  iKir  +  > 
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Table  2.1 


Guidelines  for  Stable 
Time  Step  in  Explicit  Calculations 


max 


Spring  + 

lumped  masses 


Square  of  Frequency  w 

max 


4k 

ml  +  m2 


Spring 

uniform  mass 


Beam 

uniform  mass 


or 


192c2I 


A  l 


Same  as  beam  with 
l  =  shortest  distance 
from  node  to  opposite  edge 
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Plate 


The  internal  work  is  similarly  defined  by 


wlnt’°  =  0 


(2.03) 


wint,j  =  wint,j-l  +  i  (*(e),J  .  ‘(O.j-lj  (£(e),j  +  g(e),j-l 

2  lA  iA  lA  lA 

Energy  balance  then  requires  that 


Wint’j  +  Tj  -  W6Xt>j  <  e(Wint’j  +  Tj) 


(2.84) 


If  e  is  greater  than  about  0.05,  a  solution  is  considered  unacceptable. 
The  source  of  the  energy  error  may  be  an  arrested  or  incipient  numerical 
instability,  or  excessive  truncation  error. 


2 . 6  Contact-Impact  Formulation 

Explicit  time  integrations  are  easily  modified  to  handle  the 
contact,  impact,  and  subsequent  release  of  two  bodies  when  a  diagonal, 
lumped  mass  matrix  is  used.  The  original  continuum  then  becomes  a 
collection  of  particles  (masses)  with  interacting  forces  (springs). 
This  collection  may  then  be  analyzed  by  assuming  that  linear  momentum 
is  conserved  during  collision  and  modifying  the  computed  accelerations 
and  velocities  appropriately. 

We  will  first  discuss  the  fundamental  assumptions  made  in  our 
impact  formulation.  This  done,  we  will  then  detail  the  flow  and  basic 
operation  of  the  subroutine  which  implements  this  formulation. 

The  following  assumptions  are  made: 

1.  Linear  momentum  is  conserved  during  impact. 

2.  The  actual  impact  is  smeared  over  the  time  step  At  and  is 
plastic. 

3.  The  time  step  At  is  small  enough  so  that  contact  may  take 
place  with  only  one  element  during  that  time. 
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4.  The  time  step  At  is  small  enough  so  that  any  wave  propa¬ 
gated  by  the  impact  only  reaches  nodes  directly  adjoining 


the  contact  elements;  this  is  guaranteed  in  an  explicit 
integration  by  the  Courant  stability  conditions. 

5.  Internal  (spring)  forces  remain  constant  during  impact. 

6.  The  time  step  At  is  small  enough  so  that  the  deformed 
configuration  before  impact  is  the  same  after  impact. 

7.  The  motion  of  an  adjoining  node  of  two  elements  contacted 
during  a  time  step  is  equivalent  to  a  weighted  average  of 
one  body  contacting  it. 

8.  Thermal  energy  of  impact  is  dissipated  completely  into  the 
environment  surrounding  the  bodies. 

9.  All  sliding  is  frictionless. 

10.  Rigid  body  rotation  of  the  pilot’s  head  during  the  time 
step  impact  is  negligible  when  compared  to  the  translational 
velocities. 

11.  The  plate  elements  are  treated  as  flat,  triangular  surfaces. 

The  first  assumption  is  a  fundamental  axiom  in  dynamics,  the 
second  is  a  smoothing  approximation  with  the  view  that  as  At  +  0  the 
response  approaches  a  suitable  solution.  The  third  through  fifth 
assumptions  permit  a  total  uncoupling  of  the  system  during  impact  and 
allow  one  to  work  with  a  single  element,  and  the  sixth  is  introduced 
to  avoid  any  iteration  in  determining  the  "true"  impact  point  (see 
Hughes  (1975)  for  a  treatment  of  this).  The  error  this  approximation 
induces  is  small  and  has  been  found  to  have  little  effect  on  the  results. 
The  seventh  assumption  is  needed  to  prevent  any  iteration  in  the  impact 
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analysis.  It  is  felt  that  such  an  iteration  would  be  very  harmful  to 
computing  efficiency.  The  eighth  is  made  since  the  constitutive  model 
of  the  canopy  does  not  incorporate  any  temperature  softening  or  temper¬ 
ature  melting  law  and  since  there  is  no  capability  for  a  heat  transfer 
analysis  at  present  in  the  program.  The  last  assumptions  are  made  for 
the  sake  of  smoothness  and  simplicity  of  programming. 

We  will  now  describe  the  algorithm  which  implements  contact-impact. 
We  define  here  as  master  nodes  the  canopy  nodes,  slave  nodes  the 
nodes  of  the  bird  or  helmet.  Impact  is  defined  the  first  time  step 
during  which  contact  between  two  bodies  is  detected. 

1 .  Update  all  nodal  displacements  without  considering  the 
possibility  of  contact. 

2.  Check  for  contact  between  different  bodies.  If  there 
is  no  penetration  or  contact,  then  continue  to  the 
next  time  step. 

3.  If  contact  between  a  bird  and  the  canopy  has  taken  place, 
determine  the  elements  contacted  by  each  node  of  the 
bird  and  interpolate  the  mass  of  the  bird  to  the  ad¬ 
joining  nodes  of  the  element.  Then  adjust  the  component 
of  the  displacement  of  the  bird  normal  to  the  canopy,  so 
that  it  lies  in  the  plane  of  the  plate  element  and 
satisfies  any  boundary  conditions.  If  contact  takes  place 
between  the  pilot's  helmet  and  the  canopy;  adjust  the 
vertical  displacement  of  the  head  so  that  only  one  point 
of  the  helmet  lies  on  the  surface  of  the  plate  element, 
and  interpolate  the  mass  of  the  helmet  to  the  nodes  of 
the  contact  element. 
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4.  In  the  case  of  impact,  use  conservation  of  momentum  and 
the  notion  of  plastic  impact  to  modify  the  normal 
velocities  of  master  nodes.  For  both  impact  and  con¬ 
tinuing  contact,  modify  the  normal  accelerations  to  reflect 
the  kinematic  coupling. 

5.  Obtain  the  new  slave  node  velocities  and  accelerations 
by  interpolation  on  the  element  surface. 


The  equations  used  are  as  follows.  To  check  for  contact  between 
a  node  T  and  an  element  e  with  nodes  P  =  1,  2,  3,  one  first  needs  to 
check  the  normal  distance  between  node  T  and  element  e.  Observing 
Figure  2.4,  we  see  that  this  amounts  to  calculating, 

B  =  (xT  -  Xp)  •  n 


(2.85) 


where  n  is  the  "inward"  normal  of  the  element  and  and  Xp  are  the 
coordinates  of  nodes  T  and  P;  B  is  the  distance  the  nodes  has  penetrated 
element  e,  so  if  B  >  0  the  node  T  may  have  penetrated  element  e. 

In  order  to  insure  penetration  has  occurred  for  element  e  instead 
of  another  element,  we  compute  the  projected  areas  using  the  vector 
cross-product 


Ap  2  '  ^XP+1 


-  xp+2)  x 


(XT  "  XP+2) 


(2.86) 


Here  P  +  1  denotes  the  node  on  the  counterclockwise  side  of  P  in  element 

3 

e.  Then  if  B  >  0  for  any  P  and  0<A  <  A  =  T  A,  contact  has  taken 

-  P  -  p=l  P 

place  during  the  time  step  with  element  e. 

If  contact  or  penetration  has  taken  place,  the  new  position  of  the 
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bird  node  is  determined  by 


(2.87) 


The  superscript  (+)  denotes  the  time  immediately  after  contact. 

A  slightly  different  procedure  for  locating  the  new  position  of 
the  contact  node  must  be  followed  when  node  T  has  to  satisfy  boundary 
conditions,  such  as  on  lines  of  symmetry.  Let  us  assume,  for  example, 
that  nodes  1  and  2  of  the  plate  element  lie  in  the  x^  -  plane  which 
is  a  line  of  symmetry.  If  we  followed  the  previous  procedure  for 

►  I  12 

calculating  x^  we  would  incorrectly  move  x^  out  of  the  x  -x  plane, 
because  the  plate  element  is  not  necessarily  normal  to  the  line  of 
symmetry.  Therefore,  we  let  n  =  (n^,  0)  and  use  Eq.  (2.87)  with 

this  normal  vector. 

We  may  use  a  similar  procedure  for  computing  the  new  position  when 
the  pilot's  helmet  and  the  canopy  impact.  The  helmet  is  a  sphere  of 
radius  R,  but  its  primary  node  is  at  its  center.  The  condition  for 
contact  or  penetration  then  becomes 


B  >  R 


(2.88) 


with  the  same  conditions  on  A  . 

P 

-H-  ,  .  ■* 

xT  =  xT  -  (B  -  R)  n 


The  point  of  contact,  x^ ,  is 


Eq.  (2.87)  is  similarly  replaced  by 

*  (2.89) 

(2.90) 


Once  contact  is  established  with  a  plate  element,  the  plate  element 
nodes  become  the  master  nodes  in  driving  the  normal  component  of  the 
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displacements.  For  this  purpose,  it  is  necessary  to  ascribe  the  nodal 
masses  of  the  slave  nodes,  i.e.  the  bird  or  head,  to  the  master  nodes. 
If  the  mass  of  node  T  is  denoted  p  the  mass  of  node  T  is  given  after 
contact  by 


Pp  =  Pp  +  pT  Ap/A 


(2.91) 


To  obtain  the  velocity  at  each  node,  we  note  that  only  the  normal 
direction  is  influenced  by  impact  in  the  absence  of  friction.  The 
normal  component  of  any  vector  is  given  by 


u  =  u  n . 
n  i  1 


(2.92) 


The  modifications  which  simulate  impact  are  then  as  follows: 


•+  PPUPn  +  PT  UTn 


(2.93a) 


PP  UPn+  PT  UTn 


(2.93b) 


The  velocities  and  accelerations  are  then  interpolated  to  the  birds 
positions  so  that  only  the  normal  direction  is  coupled  we  need.  For 
this  purpose,  only  the  linear  part  of  the  shape  functions  are  used, 
which  gives 


•+  1  f  .  •+ 

UTn  =  A  ^  UPn 


(2.94) 


••+  I  r  .  "T 

UTn  ”  A  T  ,  \  UPn 
1=1 


The  tangential  velocities  are  uneffected  by  impact  since  friction  is 
neglected. 
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In  time  steps  subsequent  to  impact  (continuing  contact),  Eqs. 
(2.93a)  and  (2.94a)  are  not  used;  instead  the  normal  velocities  are 
simply  updated  using  the  integration  formulas,  Eqs.  (2.52). 


SECTION  3 


PILOT  MODEL  AND  INJURY  CRITERION 


3 . 1  Spine  Model 

For  the  purpose  of  representing  the  pilot,  the  simplified  spinal 
model  (SSM)  described  by  Belytschko  and  Privitzer  (AMRL/TR-78-7)  is 
used.  The  more  complex  models  described  in  that  report,  such  as  the 
complex  spine  model  (CSM)  and  the  isolated  ligamentous  spine  with 
viscera  (ILSV),  which  were  developed  in  Belytschko  et  al  (1976)  were 
not  used  because: 

1.  Details  of  lower  spine  response  are  not  of  interest, 
so  any  model  which  represents  the  overall  flexural 
and  axial  response  of  the  spine  is  adequate. 

2.  The  computational  cost  and  computer  storage  required 
by  the  SSM  are  considerably  less  than  for  the  other 
two  models,  which  proved  of  advantage  since  a 
considerable  number  of  degrees  of  freedom  are  required 
to  model  a  canopy. 

In  the  original  SSM,  the  thoracolumbar  spine  is  modelled  by  three 
beam  elements,  T1-T10,  T10-L3,  and  L3-S1,  as  shown  in  Fig.  3.1.  The 
longitudinal  actions  of  the  viscera  are  modelled  with  four  spring 
elements  and  three  masses  between  the  pelvis  and  T10.  The  interaction 
between  the  viscera  and  the  spine  is  modelled  with  a  spring  element 
at  each  visceral  mass.  However,  the  neck  portion  of  the  model,  which 
was  modelled  in  the  original  SSM  model  by  a  single  beam  element,  has 
been  refined  by  subdividing  the  single  beam  into  three  beam  elements; 
this  allows  the  mass  of  the  neck  to  be  more  uniformly  distributed 
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See  Fig.  3.3  for 

35,44,50  details  of  Head 


Primary  Secondary 

Element  numbers 
are  underlined 


Nodes  and  Elements  in  Simplified  Spine  Model  (SSM) 


between  C7  and  Tl.  Refinement  of  the  neck  is  considered  of  importance 


in  this  study  since  its  response  plays  a  significant  role  in  head  impact. 

The  nodal  mass  distribution  and  stiffness  properties  of  the 
elements  of  SSM  are  given  in  Tables  3.1  and  3.2,  respectively.  Both 
the  masses  and  stiffnesses  were  obtained  from  the  more  complex  models 
by  combining  the  masses  and  using  appropriate  series  formulas  for  the 
stiffnesses.  The  stiffnesses  of  the  individual  vertebral  motion  segments 
were  originally  developed  by  Schultz,  et  al  (1973'),  while  the  visceral 
stiffnesses  were  obtained  by  Belytschko  and  Privitzer  (1978)  by 
matching  impedance  data.  The  inertial  data  is  based  on  Liu  and 
Wickstrom  (1973). 

Although  the  simplified  spine  model,  unlike  the  isolated  ligamentous 
spine  model  or  the  complex  spine  model,  does  not  represent  the  indi¬ 
vidual  vertebrae,  it  was  felt  desirable  to  have  the  ALFA  PLOT  graphical 
display  procedures  developed  for  the  other  spine  models  also  available 
for  the  simplified  spine  model.  This  was  accomplished  by  developing  a 
short  post-processing  routine  for  the  program  SAM,  which  linearly 
interpolates  the  unit  vectors,  b^  and  the  positions  u^  for  each  of 
the  vertebral  bodies  which  does  not  correspond  to  a  node  in  the  SSM. 

An  illustration  of  the  resulting  graphical  display  may  be  found  in 
Figure  3.2. 


III. 2  Head  and  Injury  Model 

Available  injury  models  may  be  classified  as  follows: 

1.  Simple  severity  indices  based  on  the  magnitude  and  duration 


of  the  acceleration  time  history. 


2.  Indices  based  on  responses  of  single  degree  of  freedom 
models  which  represent  the  skull-brain  impedance. 
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Figure  3. 


2  Computer  Graphics  Display  of  Simplified  Spine  Model  (SSM) 
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Table  3.1 


Stiffness  Data  for 
Simplified  Spine  Model 


Stiffnesses 


Element 

Approximate 

Axial 

Bend ing 

U  1. 

1  &  2 

Torsional 

number 

Physical 

Reg ion 

ka 

dynes 

1  U  f  1  U 

kl  &  k2 
dyne-cm 

kfc 

dyne-cm 

2 

S  L3-S1 

5 . 68E  +  8 

3.43E  + 

8 

3.5E  +  8 

3 

S  T10-L3 

3.47E  +  8 

2.15E  + 

8 

1.87E  +  8 

4 

S  T1-T10 

1.54E  +  8 

6.67E  + 

7 

4.86E  +  7 

5 

S  C6-C7 

4.32E  +  8 

2.07E  + 

9 

1.55E  +  8 

6 

S  C3-C5 

ft 

It 

II 

7 

S  C2-C4 

M 

ft 

It 

20 

Viscera 

5.4E  +  6 

21 

II 

8 . 26E  +  6 

22 

If 

1.03E  +  7 

23 

It 

8.75E  +  6 

24 

Viscera 

5.0E  +  6 

25 

< 

Spine 

5.0E  +  6 

. 

27 

Intercon¬ 

1.0E  +  7 

28 

nectors 

1.0E  +  7 

*-  - 
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Table  3.2 


Mass  Data  for  Simplified  Spine  Model 


Approximate 

Primary  Node  Physical 

Number  I  Region 


Mass 

gm 


II 


Moments  of 
gm  -  cm1 
I 


inertia 

2 


21 


31 


28 

Seat 

- 

- 

- 

- 

29 

Pelvis 

2.940E 

+ 

4 

1.308E 

+ 

6 

2.029E 

+ 

6 

1.938E 

+ 

30 

S  L3-L5 

1.660E 

+ 

3 

3.996E 

+ 

4 

5.786E 

+ 

4 

2.800E 

+ 

31 

S+V  T6-L2 

6.846E 

+ 

3 

5.328E 

+ 

5 

8.560E 

+ 

5 

8.307E 

+ 

32 

T  C6-T5 

5.896E 

+ 

3 

2.835E 

+ 

5 

6.625E 

+ 

5 

7.429E 

+ 

33 

T  C4-C5 

6.654E 

+ 

2 

1.428E 

+ 

4 

2.410E 

+ 

4 

3.861E 

+ 

34 

T  C1-C3 

6.374E 

+ 

2 

1.421E 

+ 

4 

1.263E 

+ 

4 

2.106E 

+ 

35 

Brain 

3.940E 

+ 

3 

1.420E 

+ 

5 

1.420E 

+ 

5 

1.420E 

+ 

36 

V  L4-L5 

4.295E 

+ 

3 

1.272E 

+ 

5 

2.985E 

+ 

5 

3.158E 

+ 

37 

V  L2-L3 

3.368E 

+ 

3 

1.203E 

+ 

5 

3.484E 

+ 

5 

3.746E 

+ 

38 

V  T10-L1 

5.134E 

+ 

3 

2.576E 

+ 

5 

5.598E 

+ 

5 

6.257E 

+ 

41 

Seatback 

- 

- 

- 

- 
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Skull 

6.570E 

+ 

2 

3.790E 

+ 

4 

3.790E 

+ 

4 

3.790E 

+ 

50 

Helmet 

1.400E 

+ 

3 

1.000E 

+ 

5 

1.000E 

+ 

5 

1.000E 

+ 

6 

4 

5 
5 
4 

4 

5 
5 
5 
5 

4 

5 


V  =  Viscera 
S  =  Spine 
T  =  Tbrso  or  neck 
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3.  Continuum  models  of  the  skull-brain  system  which  are  solved 
either  by  analytic  or  finite  element  methods. 

We  will  review  the  applicability  and  appropriateness  of  various 
models  for  this  study,  beginning  with  the  last  category.  Analytical 
continuum  models  of  the  head  have  been  quite  numerous  and  are  reviewed 
by  Y.  K.  Liu  (1978)  and  by  Engin  and  Engin  (1975).  Although  these 
models  have  brought  about  some  understanding  of  the  phenomological  aspects 
of  head  injury  in  impact  environments,  they  have  not  been  successfully 
correlated  with  injury.  Recent  finite  element  models  of  the  head  and 
brain  seem  to  be  substantially  more  successful,  in  that  correlations  have 
been  obtained  between  injury  and  predicted  maximum  intercranial  pressure 
(Ward  (1978)).  However,  these  models  are  so  complex  that  their 
incorporation  in  a  bird  simulation,  where  the  canopy,  head  and  spine 
must  be  simultaneously  modeled,  would  entail  exorbitant  computational 
resources . 

For  these  reasons,  our  attention  was  restricted  to  injury  models 
in  the  first  and  second  categories.  In  the  first  category,  the  most 
well  known  are  the  Wayne  State  Tolerance  Curve  (WSTC),  Lissner  (1960), 
the  Gadd  Severity  Index  (GSI),  and  the  Head  Injury  Criterion  (HIC) 
proposed  by  J.  Versace  (1971)  and  then  modified  by  NHTSA.  In  the  WSTC, 
the  injury  potential  is  based  simply  on  the  average  acceleration  and  its 
duration.  A  curve  relating  the  threshold  average  acceleration  to  the 
pulse  duration  is  given;  below  this  curve,  at  most  cerebral  concussions 
without  aftereffects  are  predicted. 

In  the  GSI,  the  injury  potential  is  expressed  in  terms  of  the 


acceleration  history  a(t)  and  duration  t  by 
*  f*  ^  2 

GSI  =  an  (t)  dt  (3.1) 

»  c 

where  n  is  a  weight  factor,  with  the  recommended  value  of  2.5.  A  GSI 
value  of  1,000  is  considered  life  threatening.  The  HIC  is  a  modifi¬ 
cation  of  the  GSI  in  which  the  following  formula  is  used 


HIC 


1  f*2 

-21^21  ^ 


a  dt 


(3.2) 


where  t^  and  t ^  ^  C2  '  C1  are  c^°8en  as  to  maximize  the  HIC.  This 
criterion  has  been  normalized  against  injury  data  so  that  an  HIC  threshold 
for  life  threatening  injury  is  1,000. 

The  second  category  of  injury  indices  are  based  on  single  degree 
of  freedom  models  of  the  skull-brain,  which  can  be  described  by  the 
equation 


z  +  20  n  z  +  o2  z  =  0  (3.3) 

where  z  represents  a  relative  displacement  or  an  effective  strain  for 
the  system,  fl  is  a  period  for  the  system  usually  based  on  impedance 
data,  and  S  the  fraction  of  critical  damping.  Several  models  of  this 
type  have  been  introduced  by  the  Vienna  Institute.  In  the  work  of  A. 
Slattenscheck  (1968),  a  natural  angular  frequency  of  634  radians 
per  second  and  a  critical  damping  coefficient  of  1.0  is  used.  Injury 
potential  is  based  on  the  maximum  value  of  z,  as  compared  to  a  toler¬ 
ance,  thus  leading  to  a  J  index  defined  by: 

J  -  (3.4) 

Zt 

where  z  equals  to  0.092  inches.  The  threshold  value  for  J  for  human 
tolerance  is  1.0.  Subsequent  models  of  this  type  have  been  based  on 
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revised  values  of  ft,  6  and  z .  The  latest  model  reported  in  an  NHSTA 
(1969)  report  uses  a  natural  frequency  of  175  radians  per  second,  a 
fraction  of  critical  damping  of  0.4  and  a  z ^  value  of  1.25  inches. 

Stalnaker,  et  al  (1971),  McElhaney,  et  al  (1973)  have  revised 
these  models  and  performed  extensive  correlations  to  injury  data.  This 
model  has  come  to  be  known  as  the  maximum  strain  criterion  (MSC)  model 
for  the  injury  potential.  Injury  is  related  to  the  displacement  of  the 
one  degree  of  freedom  system,  which  is  conceptualized  as  an  average  strain 
in  the  brain. 

In  the  MSC,  experimentally  determined  head  injury  data  for  living 
primates  is  used  to  define  a  tolerable  mean  strain  for  humans. 

McElhaney  et  al  defined  the  mean  strain  as  the  displacement  of  one  side 
of  the  head  relative  to  the  other,  divided  by  the  distance  across  the 
cranium.  Using  experimental  data  they  report  that  their  preliminary 
mean  strain  level  of  0.0061  may  be  used  as  an  injury  criterion  in 
modeling  human  head  impact,  and  called  this  a  tolerable  mean  strain, 
which  is  the  mean  "strain  corresponding  to  accelerations  required  to 
produce  an  ESI  (Estimated  Severity  of  Injury)  of  level  class  3, moderate 
but  reversible  type  of  injury." 

Put  simply,  the  MSC  states  that  an  element  strain  (for  the  elements 
representing  the  interaction  of  the  brain  and  the  skull)  less  than  0.0061 
indicates  at  most  moderate  but  reversible  injury  levels,  while  an 
element  strain  greater  than  0.0061  indicates  non-reversible  and  possibly 
fatal  injury  levels.  A  strain  value  of  0.00329  was  established  by 
Stalnaker,  et  al  as  a  survivable  acceleration  pulse. 

McElhaney,  et  al  compared  their  MSC  injury  criterion  with  other 
head  impact  injury  criteria  -  the  GSI  (Gadd  Severity  Index) ,  the  HIC 


(Head  Injury  Criterion),  the  JTI  (Vienna  Institute  Index),  the  EDI 
(Effective  Displacement  Index),  and  the  RBM  (Revised  Brain  Model).  They 
give  brief  descriptions  of  these  criteria,  and  present  the  results  of  a 
study  which  demonstrates  that  the  MSC,  EDI,  RBM  and  JTI  predict  essential¬ 
ly  the  same  injury  levels,  while  the  HIC  and  GSI  predict  somewhat  higher 
injury  levels. 

Because  of  the  computational  simplicity  and  extensive  correlation 

with  injury  data  available  for  the  MSC,  it  was  chosen  for  use  in  this 

model.  We  will  now  describe  the  way  in  which  this  head  injury  model 

was  incorporated  in  the  biodynamic  model. 

Figure  3.3  depicts  the  lateral  and  top  views  of  the  head  model. 

This  model  consists  of  three  concentric  spheres,  which  are  interconnected 

by  springs,  and  represent  the  brain,  skull  and  helmet,  respectively. 

This  model  was  incorporated  into  the  Simplified  Spine  Model  (SSM). 

Although  it  is  a  rather  simple  head  model,  it  allows  for  the  direct 

application  of  the  Mean  Strain  Criterion  (MSC). 

The  spring  constants  and  masses  for  this  model  are  based  partially 

on  mechanical  impedance  results  for  the  head  presented  by  McElhaney, 

et  al.  The  dominant  peak  for  both  lateral  and  anterior-posterior 

3 

motion  was  reported  to  be  in  the  vicinity  of  300  Hz  (1.88  x  10  rad/sec) . 

It  is  assumed  that  this  is  the  natural  frequency  of  the  brain/skull 
system.  The  natural  circular  frequency  of  a  two  mass-spring  system  is 
given  by 

,  (m.  +  m„)  2k 

ST  =  — - - - -  (3.5) 

ml  m2 

where  m^  and  m2  are  the  masses  and  k  is  an  effective  spring  constant. 


Head  with  helmet  and  brain  for  MSC  injury  model. 


I 


Solving  for  k  we  have 


m2  £2 


2  (m^  +  m2) 


McElhaney,  et  al  report  the  following  values:  skull  mass,  m^  =  2.72  x  10 

3 

grams,  brain  mass,  =  4.54  x  10  grams,  skull-brain  stiffness,  2k  = 

9 

3.45  x  10  dyne/cm.  Substituting  these  values  into  Eq.  (3.5)  yields  a 

3 

circular  frequency  of  3.67  x  10  rad/sec  (580  Hz).  This  is  approximately 
twice  the  frequency  obtained  from  the  impedance  curves.  Rather  than 
using  these  values,  we  decided  to  maintain  the  300  Hz  value,  approximate 
m^  and  m2  from  the  other  data  and  then  compute  a  spring  constant  k  from 
Eq.  (3.6). 

The  other  values  are  based  on  Hodgson  and  Patrick  (1973),  who 

developed  a  head  model  for  impact  studies.  They  report  average  values  of 
3 

4.6  x  10  grams  for  the  total  head  weight,  and  19  and  16  cm  for  the 
anterior-posterior  and  lateral  diameters,  respectively.  Hodgson  also  shows 
a  sagittal  cross-section  of  the  human  skull  which  he  used  as  his  model 
pattern.  The  mean  skull  thickness  appears  to  be  on  the  order  of  0.5  cm. 
From  this  information,  the  skull  mass,  m^,  and  brain  mass,  m2,  were 
computed  as  follows.  If  we  treat  the  brain  and  skull  as  concentric  spheres 
with  radii  of  9.5  and  10  cm  respectively,  then  the  volumes  V  are  given  by 

V  (brain)  =  y  it  (9. 5) 3 


V  (skull)  =  \  it  [103  -  9.53] 


V  (total)  =  j  it  ( 10) 3 
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Hence  the  ratios  are 


V  (brain) /V  (total)  -  0.857 

V  (skull)/V  (total)  =  0.143 

Assuming  a  constant  density  for  the  skull  and  the  brain  and  a  total 

3 

head  mass  of  4.6  x  10  grams,  and  m 2  are  given  by 

=  (.143)  (4.6  x  103)  =  6.57  x  102  grams 

m^  =  (.857)  (4.6  x  103)  =  3.94  x  103  grams 

q 

With  these  values  for  and  m,,,  Eq.  (3.6)  yields  a  value  of  1.00  x  10 
dyne/cm  for  the  springs  connecting  the  two  masses. 

The  rotational  inertias  of  the  brain  and  skull  were  computed  as 
follows.  The  rotational  inertias  of  a  solid  sphere,  I  ,  and  a  hollow 

D 

sphere,  I  ,  are  given  by 
H 

Ig  =  2mr2/5  (3.8a) 

I„  =  8mir  (r  5  -  r,5)  /  15  V  (3.8b) 

ri  0  1 

Eqs.  (3.8a)  and  (3.8b)  yield  values  of  1.42  x  103  and  3.79  x  10^  gram-cm2 
for  the  brain  and  skull  respectively. 

The  values  used  for  the  translational  and  rotational  masses  of 
the  helmet  were  those  reported  by  Belytschko,  et  al  (1976), 

3 

m  =  1.40  x  10  grains 

5  2 

I  =  1.00  x  10  gram-cm 

The  spring  constants  for  the  springs  connecting  the  helmet  to  the 
head  were  determined  as  follows.  The  material  of  the  interior  of  the 
helmet  was  taken  to  be  styrofoam.  A  specific  value  for  the  Young's 


Table  3.3;  Stiffness  and  Intertial  Properties  of  Brain/Skull/Helmet 


Head  Impact  Model 


Linear 
Stiffnesses 
dyne /cm 

Translational 
mass  (grams) 

Rotational 
mass  (gram 

Brain 

1.00  x  109 

3.94 

x  103 

1.42 

x  105 

Skull 

rigid 

6.57 

x  102 

3.79 

4 

x  10 

He lme  t 

o(D 

4.48  x  10 

1.40 

x  103 

1.00 

x  105 

(1) 


stiffness 


in  compression  only 
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modulus  is  not  available  because  styrofoam  is  a  porous  styrene  matrix 
and  its  mechanical  properties  are  functions  of  density. 

A  cubic  load-deflection  relationship  was  used 


f  =  (k  +  K  102  62)  6 


(3.9) 


The  Young’s  moduli  for  synthetic  rubbers  have  a  range  of  approximately 

500-10,000  psi.  Based  on  this,  the  Young's  modulus  was  taken  to  be  500  psi 
7  2 

(1.75  x  10  dyne/cm  ).  The  linear  stiffness  of  each  of  the  6  helmet 
spring  elements  is  then  given  by 

AF  9 

k  =  —  =  4.48  x  10  dyne/cm 

oL 

where  A  is  total  surface  area  of  the  middle  surface  of  the  helmet  which 
3  2 

is  1.52  x  10  cm  .  Table  3.3  summarizes  the  stiffness  and  inertial 
properties  of  the  brain-skull-helmet  (BSH)  system. 

From  one  dimensional  modal  analyses  of  the  SSM  with  the  brain-skull- 
helmet  (BSH)  head  impact  model,  it  was  determined  that  the  natural  frequen¬ 
cy  of  the  brain-skull  representation  is  approximately  300  Hz  for  motion 
in  all  three  global  coordinate  directions.  The  stiffness  proportional 
damping  parameter  for  the  brain  elements  was  based  on  50%  damping  of 
this  frequency, 


2  n  (300) 


5.31  x  10 


The  damping  ratio  reported  by  McElhaney,  et  al  was  approximately  0.20, 

while  the  damping  ratios  used  by  some  other  investigators  were  0.40, 

0.707  and  1.0  (as  reported  by  McElhaney,  et  al) . 

In  order  to  study  the  behavior  of  this  model,  several  simulations 

were  performed  on  the  head-spine  model  alone  with  initial  velocities 


57 


prescribed  to  the  helmet.  For  the  first  simulation,  an  initial  velocity 
vector  of 

Vo  =  (-0.707  i  -0.707  k)  (4.47  x  103)  cm/sec 

(|Vq|  =  100  mph)  was  prescribed  at  the  helmet;  the  head  was  2  in.  (5.08  cm) 
forward  from  the  seatback;  the  SSM  was  fixed  at  the  pelvis;  a  restraint 
system  consisting  of  two  seatbelts  between  the  seatback  and  the  pelvis 
was  included;  the  solution  time  increment  was  2.0  x  10  3  sec  and  600 
time  steps  were  used  for  a  duration  of  0.012  sec. 

Tables  3.4  and  3.5  list  the  peak  values  of  some  of  the  motion  and 
force  responses  respectively.  For  comparison.  Table  3.6  lists  some  of 
the  peak  motion  and  force  responses  from  a  similar  simulation  with  the 
SSM  with  a  single  rigid  body  representing  the  head  and  helmet.  Figure  3.4 
depicts  the  SSM  configurations  at  0,  4,  8,  and  12  milliseconds. 

A  comparison  of  the  results  given  in  Tables  3.4  and  3.5  demonstrate 
that  considerably  more  energy  is  transmitted  to  the  body  through  the 
rigid  head/helmet  than  through  the  BSH.  This  is  to  be  expected  since 
some  of  the  impact  energy  goes  into  deformation  of  the  BSH  elements. 

The  remaining  five  impact  simulations  consisted  of  varying  the 

initial  helmet  velocities  through  magnitudes  of  2235,  1118,  559,  279 

and  140  cm/sec  (50,  25,  12.5,  6.25  and  3.125  mph  respectively).  These 

9  9  8 

initial  helmet  velocities  correspond  to  2.62  x  10  ,  1.24  x  10  ,  6.11  x  10  , 
3.06  x  108  and  1.53  x  108  dyne  forces,  respectively  (5890,  2788,  1374,  688 
and  344  lbs,  respectively.)  These  were  short  duration  (0.0002  sec.) 
simulations  since  the  quantities  of  primary  interest  were  the  peak  strains 
in  the  springs  representing  the  skull-brain  system,  which  occur  early 
in  the  simulation.  Table  3.6  summarizes  the  brain  element  strains  for  all 


the  simulations.  The  numbers  given  are  the  averages  of  the  peak  strains 
in  the  top,  front,  bottom  and  rear  brain  elements. 

The  strains  vary  almost  linearly  with  initial  helmet  velocity  or 
helmet  impulse,  particulary  for  lower  levels.  Linear  interpolation 
between  the  lowest  two  values  results  in  an  initial  helmet  velocity  of 

g 

210  cm/sec  (4.70  mph)  and  a  helmet  force  of  2.74  x  10  dyne  (617  lbs) 
for  a  strain  equal  to  the  MSC  value  of  0.0061.  These  two  values  represent 
the  cut-off  levels  between  reversible  and  irreversible  brain  damage  as 
predicted  by  the  present  BSH  model.  It  should  be  noted  that  these  values 
are  most  likely  conservative  since  the  elasticity  of  the  brain  was 
represented  by  linear  springs.  In  reality,  the  stiffness  of  the  brain 
tissue  would  increase  with  increasing  deformation  resulting  in  smaller 
strains  that  those  predicted  by  a  linear  spring  model. 
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Table  3.4:  Peak  Motion  Responses  for  Brain-Skull-Helmet  Impact  Simulation, 
Vq  =  (.707  j  -  .707  k)  4.47  x  10^  cm/sec 


Peak  Motion  Response 


Helmet: 


y 

displacement 

5.80 

cm 

z 

displacement 

-1.37 

cm 

X 

rotation 

0. 66 

rad 

y 

acceleration 

-8.68 

X 

106 

cm/sec7 

z 

acceleration 

9.57 

X 

106 

cm/sec7 

X 

rotational  acceleration 

7.02 

X 

104 

rad/sec 

Skull: 

y 

displacement 

5.80 

z 

displacement 

-4.04 

X 

rotation 

-0.35 

y 

acceleration 

1.47 

X 

107 

z 

acceleration 

-1.36 

X 

io7 

X 

rotational  acceleration 

-2.18 

X 

105 

Brain: 

y 

displacement 

5.87 

z 

displacement 

-4.09 

X 

rotation 

0.0 

y 

acceleration 

7.27 

X 

io5 

z 

acceleration 

-6.31 

X 

io5 

X 

rotational  acceleration 

0.0 
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Table  3.4:  Continued 


Peak  Force  Response 


T1-C6 : 


axial  force 

-1.23 

X 

9 

10  dyne 

shear  force 

-2.60 

X 

10®  dyne 

8 

sagittal  moment 

7.80 

X 

10  dyne -cm 

C6-C4 : 

axial  force 

-1.42 

X 

io9 

shear  force 

-3.82 

X 

108 

sagittal  moment 

-1.20 

X 

io9 

C4-Skull: 

axial  force 

-2.62 

X 

109 

shear  force 

1.65 

X 

io9 

sagittal  moment 

-4.18 

X 

io9 

Helmet-Skull  element  forces  (dyne) 

top 

-1.59 

X 

io10 

front 

-1.48 

X 

io10 

bottom 

-2.36 

X 

io9 

rear 

-2.62 

X 

io9 

Skull -Brain  element  forces  (dyne) 

top 

-1.20 

X 

io9 

front 

-1.40 

X 

io9 

bottom 

1.28 

X 

io9 

rear 

1.51 

X 

io9 
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Table  3. A:  Continued 


cm 

Skull-Brain  element  strains  (  /cm) 

top 

front 

bottom 

rear 


-.0748 

-.1058 

.0861 

.1122 
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Table  3.5:  Peak  Force  And  Motion  Responses  For  Rigid  Head  Impact 


Simulation, 

VQ  =  (.707  j  -  .707  k)  4.47  x  10^  cm/sec 


Head/Helmet : 

y  displacement 
z  displacement 
x  rotation 
y  acceleration 
z  acceleration 
x  rotational  acceleration 


6.84  cm 
-16.60  cm 
-.42  rad 

6  2 

-4.16  x  10  cm/sec 
9.89  x  10^  cm/sec^ 
-3.24  x  104  rad/sec^ 


T1-C6: 

axial  force 
shear  force 
sagittal  moment 

C6-C4: 

axial  force 
shear  force 
sagittal  moment 


9 

-2.08  x  10  dyne 

-7.13  x  108  dyne 
9 

2.80  x  10  dyne-cm 

-2.27  x  109 
-8.73  x  108 
2.08  x  109 


C4-Head: 

axial  force 
shear  force 
satittal  moment 


-4.36  x  109 
3.13  x  109 
-8.20  x  109 
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Figure  3.4  Initial,  0.004,  O.OOS.^and  0.01^  sec  SSM  configurations 
initial  helmet  velocity  -  (-0.707  i  -0.707  k)  4.47  x  103  cm/sec 
(100  mps);  helmet  impulse  -  6.44  x  109  dynes  (1.45  x  10^  lbs). 


Table  3. 6:  Strains  Occurring  In  Brain  Elements 


Initial  Helmet  Velocity  Helmet  Impulse  Average  Peak  Brain 

F  Element  Strain 

o 


cm/ sec 

mph 

dyne- 

•sec 

lbs-sec 

4470 

100 

6.44 

X 

io9 

1.45  x  104 

0.0947 

2235 

50 

2.62 

X 

109 

5890 

0.0506 

1118 

25 

1.24 

X 

109 

2788 

0.0264 

559 

12.5 

6.11 

X 

108 

1374 

0.0135 

279 

6.25 

3.06 

X 

108 

688 

0.0068 

140 

3.125 

1.53 

X 

108 

344 

0.0034 
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SECTION  4 


CANOPY  MODEL 


4. 1  Elastic  Plastic  Material  Law 

The  initial  specifications  provided  by  WPAFB  defined  the  material  as 
a  Ramberg-Osgood  elastic-plastic  material.  Since  an  elastic-plastic  material 
with  a  piecewise  linear  hardening  curve  was  already  available  in  SADCAT, 
the  Ramberg-Osgood  curve  was  approximated  by  straight  line  segments,  so 
making  possible  the  use  of  the  available  stress-strain  law. 

In  this  stress  strain  law,  plastic  yielding  is  defined  by  a  Mises 
yield  condition  F,  which  in  plane  stress  is 


where  E  is  Young's  modulus  and  v  is  Poisson's  ratio.  The  plastic  strains 
are  given  by 
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(4.4) 


(epj  -  \  H  (a  )  {g} 


0  if  F  S  or  F  <  0  (elastic) 
2 

1  if  F  ■  J  and  F  ^  0  (plastic) 


f  ■) 

6F  /&7X 

r  -n 

2a  -  a 
x  y 

bF/te 

y 

>  -  < 

2a  -  a 
y  x 

bF/te 

xy 

V—  — 

6a 

^  xy  J 

(4.5) 


and  H  is  the  hardening  function;  by  prescribing  H  to  be  a  function  of 
only  the  yield  stress,  we  restrict  ourselves  to  isotropic  hardening. 

The  incremental  relation  between  stress  and  strain  is  then  given  by 

*  [cP]  -jAej  if  plastic  (4.6) 

■[Aoj'  »  [c]  -[as}  if  elastic 


where 


[c]  {g}t{g}  [C] 


H 


[cp]  -  [C] 

H  -  {g}T[C]  {g;  +  (2f )3  E 


(4.7) 

(4.8) 


where  is  the  current  plastic  modulus  for  the  uniaxial  stress-strain 


curve. 


4.2  Nonlinear  Elastic  Law 

An  elastic-plastic  material ,when  subjected  to  large  displacements, 
maintains  permanent  large  deformations  because  the  elastic  recovery 
is  very  small.  An  examination  of  Che  canopy  response  observed  in 
experiments  revealed  that  the  canopy  exhibited  little  permanent  deform¬ 
ation  even  when  the  maximum  displacements  were  very  large.  Therefore, 
a  nonlinear  elastic  law  was  incorporated  in  the  program  for  the 
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later  computations. 


In  this  law,  the  elastic  law 


W  -  [Cl  {*} 


(4.91 


is  used  when  the  deviatoric  stress  invariant,  Eq.  (4.1),  is  less  than 
the  yield  stress  a When  a ^  is  exceeded,  the  additional  stress  is 
related  to  any  additional  strain  through  a  matrix  [c  ]  given  by 


[C  ]  =■  f  [C] 


(4.10) 


where  E^  is  the  softer  modulus.  For  these  purposes,  a  portion  of  the 
stresses  governed  by  Eq.  (4.9),  0/  must  be  computed  by 


(4. Ill 


The  total  stress  is  then  given  by 


W =  Cc]  ({*}  -  kl) +  kl 


(4.12) 


where 


k}  ■ 


(4.13) 


m'1  w 


(4.14) 


where 


[cl*1  =  f 


0  2(l+u) 


(4.15) 


so,  whenever  is  exceeded 


(aJ  “  '  ~E~  +  aj  {ctC} 


(4.16) 


The  computational  procedure  for  evaluating  the  stresses  is  then: 

1.  given  a  state  of  strain,  the  trial  stress  \ aJ-  is  computed 


by  Eq.  (4.9). 
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2.  the  deviatoric  stress  invariant  associated  with  these 
trial  stresses  are  computed  by  Eq.  (4.1). 

3.  if  the  trial  stresses  do  not  exceed  ct^,  they  are  the 
correct  stresses;  otherwise  the  factor  Of  is  computed  by 
Eq.  (4.11). 

4.  the  stresses  and  strains  associated  with  the  initial 
modulus  are  computed  by  Eqs.  (4.13)  and  (4.14). 

5.  the  total  stresses  are  computed  by  Eq.  (4.12). 

4 . 3  Material  Constants  for  Canopy 

The  uniaxial  stress  strain  law  provided  by  AMRL  is  shown  in 
Fig.  4.1.  Since  the  software  used  in  this  investigation  assumed  a 
piecewise  linear  approximation  to  the  stress  strain  law,  it  was  fit 
using  the  following  constants: 

10  2 

E  ■  2.413  x  10  dynes/cm 
u  -  0.31 

8  2 

Ep  ■  3.32  x  10  dynes /cm 

8  2 

Oy  **  6.76  x  10  dynes/cm 

A  plot  of  the  piecewise  linear  approximation  to  the  uniaxial  stress- 
strain  curve  is  shown  in  Fig.  4.1  also. 

The  use  of  this  approximation  can  be  justified  by  the  fact  that  in 
simulations  where  the  canopy  contacts  the  head,  strains  are  on  the 
order  of  17.  to  57..  For  strains  of  this  magnitude,  the  small  discrepancies 
between  the  piecewise  linear  approximation  and  the  Ramberg-Osgood  fit 
should  be  of  little  importance. 
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SECTION  V 


RESULTS  OF  SIMULATIONS 


5. 1  Bird-Canopy  Simulations 

For  the  purpose  of  determining  the  mesh  refinement  needed  in  the 
canopy  and  the  validity  of  the  bird-impact  model,  we  made  a  series  of 
preliminary  analyses  with  only  the  canopy,  and  no  pilot.  In  these  sim¬ 
ulations,  the  bird  was  represented  by  a  single  pentahedral  element  with 
the  following  material  constants 

10  2 

Bulk  modulus  =  2.6  x  10  dynes/cm 

3 

Density  =  1  gm/cm 

Shear  modulus  =  0.0 

A  shear  modulus  could  not  be  included  because  the  pentahedral  element  does 
not  account  for  shear  stresses.  The  pentahedral  element  for  the  bird  was 
then  sized  so  that  the  desired  mass  was  represented. 

In  these  simulations,  the  canopy  material  was  an  uncoated  monolithic 
polycarbonate,  1.27  cm  (%  in)  in  thickness.  Strain  rate  effects  were  not 
considered  since  they  were  incorporated  in  the  material  law  provided,  which 
is  shown  in  Fig.  4.1.  The  elastic-plastic  law  was  used.  The  initial  fine 
canopy  model  is  shown  in  Fig.  5.1,  and  the  deformation  is  shown  in  Figs.  5.2 
and  5.3  for  the  coarse  and  fine  models.  In  this  case,  the  bird  mass  was 
0.98  kg  (2.16  lbs)  and  its  velocity  was  177.7  m/sec  in  the  horizontal  direction. 
The  point  of  impact  was  on  the  centerline  of  the  canopy;  location  in  the  x-z 
plane  is  shown  in  the  initial  views  of  each  of  the  canopies.  Since  the 
resulting  problem  was  symmetric,  only  half  of  the  canopy  was  modelled. 

The  maximum  deflection  for  the  3  meshes  are  given  in  Table  5.1  and 
compared  to  the  experimental  result  of  Run  GDI. 
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Figure  5.1  Top  and  Side  View  of  Undeforaed  Fine  Mesh  Model 
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Figure  5.2  Deformed  Configurations  of  Fine  Mesh  Model 
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Deformed  Configurations  of  Coarse  Mesh  Model 


a 


Table  5.1 


Comparison  of  Maximum  Deflections 


Coarse  Mixed  Fine  Experiment 

_ Mesh _ Mesh _ Mesh _ GDI _ 

Maximum  6.6  cm  7.1  cm  7.4  cm  12.7  cm 

deflection  (2.8  in)  (2.8  in)  (2.9  in)  (5.0  in) 


The  maximum  deflections  underestimate  the  experimental  value  signifi¬ 
cantly,  although  the  disparity  decreases  with  mesh  refinement.  Unfortunately, 
the  coarser  meshes  do  not  exhibit  the  wave  pattern  traveling  along  the  canopy 
and  even  the  fine  mesh  only  exhibits  the  beginning  of  a  wave  which  dies 
rapidly  as  it  progresses.  It  could  not  be  established  whether  this  disagree¬ 
ment  between  the  simulation  and  experiment  was  a  result  of  an  inadequate 
contact-impact  model  (i.e.  neglecting  friction  and  shear  in  the  bird)  or 
whether  the  canopy  models  were  not  sufficiently  refined.  The  latter  possi¬ 
bility  could  not  be  explored  because  the  refined  mesh  used  all  available 
core;  even  it  could  not  fit  with  the  pilot  model  simultaneously. 

Part  of  the  problem  was  found  to  be  due  to  the  coarseness  of  the  bird 
model.  To  illustrate  this,  we  impacted  a  single  mass  point  of  the  same  mass 
as  the  bird  (981  gm);  at  26,822  cm/sec  (521  Knots)  against  the  canopy  and 
compared  an  identical  bird  simulation  with  a  pentahedron.  The  single  point 
mass  model  of  the  bird  yielded  a  maximum  deflection  of  20  cm,  whereas  the 
pentahedral  model  yielded  a  deflection  of  12.9  cm.  Evidently,  the  momentum 
of  the  back  nodes  of  the  pentahedron  are  not  transferred  to  the  canopy 
adequately.  Furthermore,  both  models  ricochet  off  the  canopy  after  causing 
maximum  deformation,  thus  eliminating  the  wave  effect.  The  deflection 
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pattern  for  the  single  node  bird  model  is  shown  in  Fig.  5.4.  The  discrepancy 
between  experiment  and  simulation  lies  in  the  failure  of  the  plate 
elements  along  the  subsequent  path  of  the  bird  to  deform  prior  to  its 
arrival,  thus  causing  a  large  slope  in  the  subsequent  path  which  causes 
the  bird  to  ricochet.  We  were  unable  to  determine  whether  this  unfor¬ 
tunate  phenomenon  was  due  to  insufficient  discretization  of  the  bird, 
canopy  or  both.  However,  in  order  to  avoid  this  difficulty,  we  simulated 
the  bird  impact  by  prescribing  an  initial  impulse  to  the  impacted  part  of 
the  canopy  in  most  of  the  simulations  and  retained  the  contact-impact  al¬ 
gorithm  only  for  modelling  canopy-helmet  impact. 

5.2  Canopy  Impulse-Pilot  Simulations 

Only  one  simulation  was  successfully  completed  with  a  complete  bird, 
canopy  and  pilot  impact,  but  because  at  the  time  of  its  completion,  computer 
generated  graphic  displays  of  the  complete  system  were  not  available,  so 
we  will  report  it  later  in  the  section. 

The  first  simulation  obtained  with  a  complete  graphical  display  was  a 

114  kg-m/sec  impulse  applied  normal  to  the  canopy  coarse  mesh 

model, 32  cm  anterior  to  the  front  of  the  pilot's  helmet.  The  impulse 

given  here  is  normal  to  the  canopy  and  corresponds  to  a  perfectly 

plastic,  frictionless  impact  with  a  bird  in  a  horizontal  path  with  a 

momentum  of  251  kg-m/sec;  the  difference  between  the  momentum  and  the 

impulse  reflects  the  angle  of  27°  between  the  surface  of  the  canopy  and 

the  horizontal  axis.  For  a  0.98  kg  bird  (2.16  lbs),  this  corresponds 

to  an  initial  velocity  of  257  m/sec  (500  knots),  and  a  kinetic  energy 
3 

of  31.8  x  10  N-m.  For  a  1.82  kg  bird,  this  corresponds  to  an  initial 

3 

velocity  of  138  m/sec  (270  knots)  and  a  kinetic  energy  of  17.3  x  10  N-m. 


TIME  2.5  E  -  0  3 


re  5.4  High  Velocity  Bird  Impact  Simulation  with  Single  Node 
Point  Model  of  Bird 
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The  latter  corresponds  closely  to  the  USAF  "goal"  tolerance  bird 
strike  of  a  1.8  kg  bird  at  350  knots. 

The  pilot  in  these  simulations  was  positioned  so  that  the  top 
surface  of  his  helmet  was  5  cm  below  the  canopy  on  the  centerline; 
however,  as  is  apparent  from  the  front  view,  because  the  element  by 
the  centerline  is  inclined  it  only  clears  the  helmet  by  1.5  cm.  The 
pilot's  head  was  5  cm  anterior  to  the  seat-back,  so  that  initial 
lateral  motions  were  not  restrained  by  the  seat-back. 

The  configurations  of  the  model  of  the  canopy  and  pilot  are  shown 
in  Fig.  5.5.  Note  that  although  the  elements  are  depicted  as  planar, 
with  discontinuities  in  the  slope,  in  the  analysis  the  canopy  is 
treated  by  cubic  displacements  which  are  almost  completely  continuously 
differentiable.  The  maximum  deflection  for  this  simulation  is  22  cm, 
which  is  actually  somewhat  higher  than  the  maximum  reported  for  the 
experiments  by  Specker,  et  al  (1978)  for  this  energy  level;  Specker, 
et  al  reported  approximately  a  20  cm  maximum  and  15  cm  maximum  deflection 
at  the  pilot  station  for  this  energy  level.  At  the  pilot  however,  the 
maximum  deflection  obtained  here  is  quite  small  (~  8  cm)  and  the 
nodes  behind  the  pilot  hardly  deflect  at  all;  this  is  a  consequence 
of  the  impulse  simulation  of  bird  strike.  Thus,  the  results  should 
be  interpreted  more  in  terms  of  the  interference  between  the  pilot's 
helmet  and  canopy  (which  is  about  7  cm).  The  MSC  results  are  given 
in  Table  5.2;  the  maximum  strains  in  the  vertical  element  and  the 
maximum  strains  in  the  anterior-posterior  element  are  given.  Lateral 
element  strains  were  always  almost  zero  because  of  symmetry.  The  MSC 
value  represents  the  maximum  of  the  square  root  of  the  sum  of  the  strains. 
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Table  5.2 


Maximum  Strains  Criterion  for  Simulations 


S imulation 
Number 

Canopy 
Impulse 
kg  m/sec 

Notes 

Maximum 

Vertical 

Strains 

Maximum 
Ant. -Post. 
Strains 

Maximum 

Strain 

(MSC) 

1 

114 

-5 . 97o 

3.07. 

6.47. 

2 

144 

reduced  canopy 
strength 

-6.07. 

3.27. 

6.77. 

3 

85.5 

-1.37. 

0.57, 

1.47. 

4 

85.5 

pilot  raised 

1.0  cm 

-3.77. 

1.07. 

4.07. 
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This  simulation  showed  significant  rotations  of  the  pilot's  head. 

This  can  be  observed  from  the  earline,  which  is  coincident  with  the 
vertical  axis  in  the  undeformed  display  and  at  0.25  m/sec.  However, 
in  the  0.63  and  1.0  m/sec  displays  the  earline  is  displaced  to  the 
back  and  rotated  counterclockwise.  Since  rotational  accelerations  are 
of  importance  in  brain  stem  injuries  and  impairment  but  are  not  reflected 
in  the  MSC,  the  application  of  this  injury  criterion  in  future  work 
on  bird-strike  should  be  re-examined. 

The  elastic-plastic  law  was  used  in  this  simulation,  so  as  can 
be  seen  from  Fig.  5.5  the  canopy  does  not  rebound  to  its  initial  shape. 

It  can  be  seen  from  Table  5.2  that  the  larger  strains  are  in  the 
vertical  direction,  reflecting  the  predominantly  vertical  direction 
of  the  impact  to  the  helmet.  The  MSC  is  above  the  MSC  criterion  for 
irreversible  injury  of  6.1%. 

Simulation  2  repeated  the  same  impulse  but  the  canopy  yield  stress 
and  moduli  that  are  given  in  Chapter  4  were  reduced  by  15%.  The 
purpose  of  this  simulation  was  to  investigate  the  sensitivity  of  the 
results  to  the  material  properties  of  the  canopy.  As  can  be  seen 
from  Table  5.2  the  change  in  the  results  is  even  smaller  than  the 
change  in  the  material  properties. 

In  simulation  3,  an  impulse  of  85.5  m/sec  was  applied  normal  to 
the  canopy  at  the  same  point  as  the  previous  simulations.  The 
position  of  the  pilot  and  the  canopy  properties  were  identical  to 
simulation  2.  However,  a  mixed  mesh  with  smaller  elements  was  placed 
in  the  vicinity  of  the  pilot,  and  the  nonlinear  elastic  law  was  used. 

This  simulation  reflects  a  bird  with  a  horizontal  momentum  of  188  kg-m/sec, 
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Figure  5.5  Response  of  Canopy  and  Pilot  for  114  kg-s/sec  Impulse  (Simulation  1) 


figure  5.5c 


Figure  5.6  Response  of  Canopy  and  Pilot  for  85.5  kg -in/ sec  Impulse  (Simulation  3) 


Figure  5 . 6b 
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Figure  5.9  Displacement  Time  History  of  Node  88  (above  anterior 
surface  of  pilot) 
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which  is  associated  with  a  velocity  and  kinetic  energy  of  192.4  m/sec 

3 

(374  knots)  and  18.1  x  10  N-m,  respectively,  for  a  0.98  kg  (2.16  lbs) 

3 

bird,  and  103.5  m/sec  (201  knots)  and  9  x  10  N-m,  respectively,  for 
a  1.82  kg  (4  lbs)  bird. 

The  graphical  display  of  this  simulation  is  shown  in  Fig.  5.6  The 
maximum  deflection  at  the  impact  point  is  8  cm;  which  agrees  reasonably 
well  with  the  reported  values  at  4  to  12  cm  for  this  energy  level  at 
the  pilot  station; the  maximum  deflection  is  about  5.5  cm,  resulting  barely 
in  interference.  In  fact,  the  impact  between  the  canopy  and  helmet  is 
not  evident  in  Fig.  5.6  because  it  occurs  off  the  centerline  in  the 
middle  of  an  element. 

The  time  histories  of  a  vertical  (top)  and  anterior  posterior 
(front)  brain-spring  element  are  given  in  Figs.  5.7  and  5.8,  respectively. 
As  can  be  seen,  the  initial  strains  are  compressive,  followed  by  a 
rebound  of  almost  the  same  magnitude;  the  rebound  appears  to  be  due 
to  the  overshoot  of  the  brain  mass,  for  it  occurs  far  before  the 
helmet  has  a  chance  to  contact  the  seatback;  in  fact,  for  the  time 
duration  of  this  simulation,  no  motion  of  the  head  is  detectable  in 
Fig.  5.6.  However,  the  maximum  strains  are  considerably  below  the 
tolerable  level. 

Figure  5.9  shows  the  displacement  time  history  of  node  88  of  the 
canopy.  It  can  be  seen  that  the  node  deflects  downward  with  a  large 
velocity  and  then  slowly  rebounds  to  its  initial  position;  this  is 
a  consequence  of  the  nonlinear  elastic  law  used  in  this  simulation. 

The  fourth  simulation  was  identical  to  the  third  but  the  pilot 
was  initially  1  cm  closer  to  the  canopy,  thus  increasing  the  expected 


interference.  The  character  of  this  simulation  was  identical  to  the 
previous  ones  but  the  MSC  was  increased  significantly.  It  can  be 
seen  from  Table  5.2  that  this  decrease  in  interference  is  sufficient 
to  raise  the  MSC  above  the  value  associated  with  no  impairment. 

5.3  Complete  Simulation  of  Bird  Strike  -  Canopy-Pilot. 

A  complete  simulation  of  bird  strike  and  pilot  response  was  made 
with  a  1.8  kg  bird  with  a  horizontal  velocity  of  250  m/sec.  The  canopy 
deformation  is  shown  in  Figure  5.10;  it  shows  the  progression  of  a 
wave-like  pattern  up  the  canopy  but  after  the  second  indentation  the 
bird  ricochets  off  the  canopy.  The  canopy  strikes  the  helmet  at 
2.75  msec,  imparting  an  initial  velocity  to  the  head  with  a  horizontal 
component  of  4.2  m/sec  and  a  vertical  component  of  -1.8  m/sec.  This 
initial  velocity  is  somewhat  different  than  those  found  in  the  impulse 
simulations,  for  its  horizontal  component  is  larger  than  the  vertical 
component.  The  simulation  was  not  continued  until  the  maximum  strain 
occurred  in  the  brain,  but  using  these  initial  conditions  it  can  be 
seen  from  the  results  of  Chapter  3  that  a  maximum  strain  of  about  1.35% 
is  expected.  This  is  substantially  below  the  tolerance  limit  without 
impairment;  the  low  value  is  probably  due  to  the  bird's  failure  to 
stay  on  the  canopy. 

The  cost  of  this  simulation  was  about  $500.00  on  the  Northwestern 
CDC6600  system,  compared  to  a  cost  of  $75.00  for  the  impulse  simulations. 
Part  of  the  additional  cost  is  due  to  the  mesh  refinement,  which 
requires  that  more  elements  be  used  and  in  addition,  that  a  substantially 
smaller  time  step  be  used.  Furthermore,  when  the  entire  bird  strike 
event  is  simulated,  the  impact  to  the  pilot  does  not  take  place  until 
much  later,  so  that  the  total  simulation  time  has  to  be  increased. 
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